home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Software Vault: The Diamond Collection
/
The Diamond Collection (Software Vault)(Digital Impact).ISO
/
cdr44
/
frasrc19.zip
/
LYAPUNOV.ASM
< prev
next >
Wrap
Assembly Source File
|
1994-04-09
|
21KB
|
511 lines
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Wesley Loewer's attempt at writing his own
; 80x87 assembler implementation of Lyapunov fractals
; based on lyapunov_cycles_in_c() in miscfrac.c
;
; Nicholas Wilt, April 1992, originally wrote an 80x87 assembler
; implementation of Lyapunov fractals, but I couldn't get his
; lyapunov_cycles() to work right with fractint.
; So I'm starting mostly from scratch with credits to Nicholas Wilt's
; code marked with NW.
.8086
.8087
.MODEL medium,c
; Increase the overflowcheck value at your own risk.
; Bigger value -> check less often -> faster run time, increased risk of overflow.
; I've had failures with 6 when using "set no87" as emulation can be less forgiving.
overflowcheck EQU 5
EXTRN Population:QWORD,Rate:QWORD
EXTRN colors:WORD, maxit:DWORD, lyaLength:WORD
EXTRN lyaRxy:WORD, LogFlag:WORD
EXTRN fpu:WORD
EXTRN filter_cycles:DWORD
.DATA
ALIGN 2
BigNum DD 100000.0
half DD 0.5 ; for rounding to integers
e_10 DQ 22026.4657948 ; e^10
e_neg10 DQ 4.53999297625e-5 ; e^(-10)
.CODE
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; BifurcLambda - not used in lyapunov.asm, but used in miscfrac.c
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
PUBLIC BifurcLambda
BifurcLambda PROC
; Population = Rate * Population * (1 - Population);
; return (fabs(Population) > BIG);
push bp ; Establish stack frame
mov bp,sp ;
sub sp,2 ; Local for 80x87 flags
fld Population ; Get population
fld st ; Population Population
fld1 ; 1 Population Population
fsubr ; 1 - Population Population
fmul ; Population * (1 - Population)
fmul Rate ; Rate * Population * (1 - Population)
fst Population ; Store in Population
fabs ; Take absolute value
fcomp BigNum ; Compare to 100000.0
fstsw [bp-2] ; Return 1 if greater than 100000.0
mov ax,[bp-2] ;
sahf ;
ja Return1 ;
xor ax,ax ;
jmp short Done ;
Return1:
mov ax,1 ;
Done: inc sp ; Deallocate locals
inc sp ;
pop bp ; Restore stack frame
ret ; Return
BifurcLambda ENDP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; OneMinusExpMacro - calculates e^x
; parameter : x in st(0)
; returns : e^x in st(0)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
OneMinusExpMacro MACRO
LOCAL not_387, positive_x, was_positive_x
LOCAL long_way_386, long_way_286
LOCAL positive_x_long, was_positive_x_long, done
; To take e^x, one can use 2^(x*log_2(e)), however on on 8087/287
; the valid domain for x is 0 <= x <= 0.5 while the 387/487 allows
; -1.0 <= x <= +1.0.
; I wrote the following 387+ specific code to take advantage of the
; improved domain in the f2xm1 command, but the 287- code seems to work
; about as fast.
cmp fpu,387 ; is fpu >= 387 ?
jl not_387
;.386 ; a 387 or better is present so we might as well use .386/7 here
;.387 ; so "fstsw ax" can be used. Note that the same can be accomplished
.286 ; with .286/7 which is recognized by more more assemblers.
.287 ; (ie: MS QuickAssembler doesn't recognize .386/7)
; |x*log_2(e)| must be <= 1.0 in order to work with 386 or greater.
; It usually is, so it's worth taking these extra steps to check.
fldl2e ; log_2(e) x
fmul ; x*log_2(e)
;begining of short_way code
fld1 ; 1 x*log_2(e)
fld st(1) ; x*log_2(e) 1 x*log_2(e)
fabs ; |x*log_2(e)| 1 x*log_2(e)
fcompp ; x*log_2(e)
fstsw ax
sahf
ja long_way_386 ; do it the long way
f2xm1 ; e^x-1=2^(x*log_2(e))-1
fchs ; 1-e^x which is what we wanted anyway
jmp done ; done here.
long_way_386:
; mostly taken from NW's code
fld st ; x x
frndint ; i x
fsub st(1),st ; i x-i
fxch ; x-i i
f2xm1 ; e^(x-i)-1 i
fld1 ; 1 e^(x-i)-1 i
fadd ; e^(x-i) i
fscale ; e^x i
fstp st(1) ; e^x
fld1 ; 1 e^x
fsubr ; 1-e^x
jmp done
not_387:
.8086
.8087
; |x*log_2(e)| must be <= 0.5 in order to work with 286 or less.
; It usually is, so it's worth taking these extra steps to check.
fldl2e ; log_2(e) x
fmul ; x*log_2(e)
;begining of short_way code
fld st ; x*log_2(e) x*log_2(e)
fabs ; |x*log_2(e)| x*log_2(e)
fcomp half ; x*log_2(e)
fstsw status_word
fwait
mov ax,status_word
sahf
ja long_way_286 ; do it the long way
; 286 or less requires x>=0, if x is negative, use e^(-x) = 1/e^x
ftst ; x
fstsw status_word
fwait
mov ax,status_word
sahf
jae positive_x
fabs ; -x
f2xm1 ; e^(-x)-1=2^(-x*log_2(e))-1
fld1 ; 1 e^(-x)-1
fadd st,st(1) ; e^(-x) e^(-x)-1
fdiv ; 1-e^x=(e^(-x)-1)/e^(-x)
jmp SHORT done ; done here.
positive_x:
f2xm1 ; e^x-1=2^(x*log_2(e))-1
fchs ; 1-e^x which is what we wanted anyway
jmp SHORT done ; done here.
long_way_286:
; mostly taken from NW's code
fld st ; x x
frndint ; i x
fsub st(1),st ; i x-i
fxch ; x-i i
; x-i could be pos or neg
; 286 or less requires x>=0, if negative, use e^(-x) = 1/e^x
ftst
fstsw status_word
fwait
mov ax,status_word
sahf
jae positive_x_long_way
fabs ; i-x
f2xm1 ; e^(i-x)-1 i
fld1 ; 1 e^(i-x)-1 i
fadd st(1),st ; 1 e^(i-x) i
fdivr ; e^(x-i) i
fscale ; e^x i
fstp st(1) ; e^x
fld1 ; 1 e^x
fsubr ; 1-e^x
jmp SHORT done ; done here.
positive_x_long_way: ; x-i
f2xm1 ; e^(x-i)-1 i
fld1 ; 1 e^(x-i)-1 i
fadd ; e^(x-i) i
fscale ; e^x i
fstp st(1) ; e^x
fld1 ; 1 e^x
fsubr ; 1-e^x
done:
ENDM ; OneMinusExpMacro
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; modeled from lyapunov_cycles_in_c() in miscfrac.c
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;int lyapunov_cycles_in_c(int filter_cycles, double a, double b) {
; int color, count, i, lnadjust;
; double lyap, total, temp;
PUBLIC lyapunov_cycles
lyapunov_cycles PROC USES si di, a:QWORD, b:QWORD
LOCAL color:WORD, i:WORD, lnadjust:WORD
LOCAL halfmax:DWORD, status_word:WORD
LOCAL total:QWORD
; for (i=0; i<filter_cycles; i++) {
; for (count=0; count<lyaLength; count++) {
; Rate = lyaRxy[count] ? a : b;
; if (curfractalspecific->orbitcalc()) {
; overflow = TRUE;
; goto jumpout;
; }
; }
; }
; Popalation is left on the stack during most of this procedure
fld Population ; Population
mov bx,1 ; using bx for overflowcheck counter
mov ax,word ptr filter_cycles
mov dx,word ptr filter_cycles+2
add ax,1
adc dx,0 ; add 1 because dec at beginning
mov si, ax ; using si for low part of i
mov i, dx
filter_cycles_top:
dec si
jne hop_skip_and_jump ; si > 0 :
cmp i,0
je filter_cycles_bottom
dec i
hop_skip_and_jump:
mov cx, lyaLength ; using cx for count
mov di,OFFSET lyaRxy ; use di as lyaRxy[] index
; no need to compare cx with lyaLength at top of loop
; since lyaLength is guaranteed to be > 0
filter_cycles_count_top:
cmp WORD PTR [di],0 ; lyaRxy[count] ?
jz filter_cycles_use_b
fld a ; if lyaRxy[count]==true, use a
jmp SHORT filter_cycles_used_a
filter_cycles_use_b:
fld b ; if lyaRxy[count]==false, use b
filter_cycles_used_a:
; leave Rate on stack for use in BifurcLambdaMacro
; BifurcLambdaMacro ; returns in flag register
fld st(1) ; Population Rate Population
fmul st,st ; Population^2 Rate Population
fsubp st(2),st ; Rate Population-Population^2
; =Population*(1-Population)
fmul ; Rate*Population*(1-Population)
dec bx ; decrement overflowcheck counter
jnz filter_cycles_skip_overflow_check ; can we skip check?
fld st ; NewPopulation NewPopulation
fabs ; Take absolute value
fcomp BigNum ; Compare to 100000.0
fstsw status_word ;
fwait
mov bx,overflowcheck ; reset overflowcheck counter
mov ax,status_word ;
sahf ; NewPopulation
ja overflowed_with_Pop
filter_cycles_skip_overflow_check:
add di,2 ; di += sizeof(*lyaRxy)
loop filter_cycles_count_top ; if --cx > 0 then loop
jmp filter_cycles_top ; let's do it again
filter_cycles_bottom:
; total = 1.0;
; lnadjust = 0;
fld1
fstp total
mov lnadjust,0
; for (i=0; i < maxit/2; i++) {
; for (count = 0; count < lyaLength; count++) {
; Rate = lyaRxy[count] ? a : b;
; if (curfractalspecific->orbitcalc()) {
; overflow = TRUE;
; goto jumpout;
; }
; temp = fabs(Rate-2.0*Rate*Population);
; if ((total *= (temp))==0) {
; overflow = TRUE;
; goto jumpout;
; }
; }
; while (total > 22026.4657948) {
; total *= 0.0000453999297625;
; lnadjust += 10;
; }
; while (total < 0.0000453999297625) {
; total *= 22026.4657948;
; lnadjust -= 10;
; }
; }
; don't forget Population is still on stack
mov ax,word ptr maxit ; calculate halfmax
mov dx,word ptr maxit+2 ; calculate halfmax
shr dx,1
rcr ax,1
mov word ptr halfmax,ax
mov word ptr halfmax+2,dx
add ax,1
adc dx,0 ; add 1 because dec at beginning
mov i,dx ; count down from halfmax
mov si,ax ; using si for low part of i
halfmax_top:
dec si
jne init_halfmax_count ; si > 0 ?
cmp i,0
je step_to_halfmax_bottom
dec i
jmp short init_halfmax_count ; yes, continue on with loop
step_to_halfmax_bottom:
jmp halfmax_bottom ; if not, end loop
init_halfmax_count:
mov cx, lyaLength ; using cx for count
mov di,OFFSET lyaRxy ; use di as lyaRxy[] index
; no need to compare cx with lyaLength at top of loop
; since lyaLength is guaranteed to be > 0
halfmax_count_top:
cmp WORD PTR [di],0 ; lyaRxy[count] ?
jz halfmax_use_b
fld a ; if lyaRxy[count]==true, use a
jmp SHORT halfmax_used_a
halfmax_use_b:
fld b ; if lyaRxy[count]==false, use b
halfmax_used_a:
; save Rate, but leave it on stack for use in BifurcLambdaMacro
fst Rate ; save for not_overflowed use
;BifurcLambdaMacro ; returns in flag register
fld st(1) ; Population Rate Population
fmul st,st ; Population^2 Rate Population
fsubp st(2),st ; Rate Population-Population^2
; =Population*(1-Population)
fmul ; Rate*Population*(1-Population)
dec bx ; decrement overflowcheck counter
jnz not_overflowed ; can we skip overflow check?
fld st ; NewPopulation NewPopulation
fabs ; Take absolute value
fcomp BigNum ; Compare to 100000.0
fstsw status_word ;
fwait
mov bx,overflowcheck ; reset overflowcheck counter
mov ax,status_word ;
sahf ;
jbe not_overflowed
; putting overflowed_with_Pop: here saves 2 long jumps in inner loops
overflowed_with_Pop:
fstp Population ; save Population and clear stack
jmp color_zero ;
not_overflowed: ; Population
fld st ; Population Population
; slightly faster _not_ to fld Rate here
fmul Rate ; Rate*Population Population
fadd st,st ; 2.0*Rate*Population Population
fsubr Rate ; Rate-2.0*Rate*Population Population
fabs ; fabs(Rate-2.0*Rate*Population) Population
fmul total ; total*fabs(Rate-2.0*Rate*Population) Population
ftst ; compare to 0
fstp total ; save the new total
fstsw status_word ; Population
fwait
mov ax,status_word
sahf
jz overflowed_with_Pop ; if total==0, then treat as overflow
add di,2 ; di += sizeof(*lyaRxy)
loop halfmax_count_top ; if --cx > 0 then loop
fld total ; total Population
too_big_top:
fcom e_10 ; total Population
fstsw status_word
fwait
mov ax,status_word
sahf
jna too_big_bottom
fmul e_neg10 ; total*e_neg10 Population
add lnadjust,10
jmp SHORT too_big_top
too_big_bottom:
too_small_top: ; total Population
fcom e_neg10 ; total Population
fstsw status_word
fwait
mov ax,status_word
sahf
jnb too_small_bottom
fmul e_10 ; total*e_10 Population
sub lnadjust,10
jmp SHORT too_small_top
too_small_bottom:
fstp total ; save as total
jmp halfmax_top ; let's do it again
halfmax_bottom:
; better make another check, just to be sure
; NewPopulation
cmp bx,overflowcheck ; was overflow just checked?
jl last_overflowcheck ; if not, better check one last time
fstp Population ; save Population and clear stack
jmp short jumpout ; skip overflowcheck
last_overflowcheck:
fst Population ; save new Population
fabs ; |NewPopulation|
fcomp BigNum ; Compare to 100000.0
fstsw status_word ;
fwait
mov ax,status_word
sahf
ja color_zero ; overflowed
jumpout:
; if (overflow || total <= 0 || (temp = log(total) + lnadjust) > 0)
; color = 0;
; else {
; if (LogFlag)
; lyap = -temp/((double) lyaLength*i);
; else
; lyap = 1 - exp(temp/((double) lyaLength*i));
; color = 1 + (int)(lyap * (colors-1));
; }
; return color;
;}
; no use testing for the overflow variable as you
; cannot get here and have overflow be true
fld total ; total
ftst ; is total <= 0 ?
fstsw status_word
fwait
mov ax,status_word
sahf
ja total_not_neg ; if not, continue
fstp st ; pop total from stack
jmp SHORT color_zero ; if so, set color to 0
total_not_neg: ; total is still on stack
fldln2 ; ln(2) total
fxch ; total ln(2)
fyl2x ; ln(total)=ln(2)*log_2(total)
fiadd lnadjust ; ln(total)+lnadjust
ftst ; compare against 0
fstsw status_word
fwait
mov ax,status_word
sahf
jbe not_positive
fstp st ; pop temp from stack
color_zero:
xor ax,ax ; return color 0
jmp thats_all
not_positive: ; temp is still on stack
fild lyaLength ; lyaLength temp
fimul halfmax ; lyaLength*halfmax temp
fdiv ; temp/(lyaLength*halfmax)
cmp LogFlag,0 ; is LogFlag set?
jz LogFlag_not_set ; if !LogFlag goto LogFlag_not_set:
fchs ; -temp/(lyaLength*halfmax)
jmp calc_color
LogFlag_not_set: ; temp/(lyaLength*halfmax)
OneMinusExpMacro ; 1-exp(temp/(lyaLength*halfmax))
calc_color:
; what is now left on the stack is what the C code calls lyap
; lyap
mov ax,colors ; colors
dec ax ; colors-1
mov i,ax ; temporary storage
fimul i ; lyap*(colors-1)
; the "half" makes ASM round like C does
fsub half ; sub 0.5 for rounding purposes
fistp i ; temporary = (int)(lyap*(colors-1))
fwait ; one moment please...
mov ax,i ; ax = temporary
inc ax ; ax = 1 + (int)(lyap * (colors-1));
thats_all:
mov color, ax
ret
lyapunov_cycles endp
END
